Control of wavepacket spreading in nonlinear finite disordered lattices 
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In the absence of nonlinearity all normal modes (NMs) of a chain with disorder are spatially local- 
ized (Anderson localization) . We study the action of nonlinearity, whose strength is ramped linearly 
in time. It leads to a spreading of a wavepacket due to interaction with and population of distant 
NMs. Eventually the nonlinearity induced frequency shifts take over, and the wavepacket becomes 
selftrapped. On finite chains a critical ramping speed is obtained, which separates delocalized final 
states from localized ones. The critical value depends on the strength of disorder and is largest 
when the localization length matches the system size. 
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I. INTRODUCTION 

Spatial discreteness and nonlinearity have been probed 
recently within the context of Bose-Einstein condensates 
(BEC) in optical lattices [ij and propagation of light 
in nonlinear optical waveguide arrays to name a 
few. The balance between these ingredients allows for 
the excitation of localized structures known as discrete 
breathers/solitons Q. Localization is due to nonlinearity 
induced frequency shifts, which tune a localized excita- 
tion out of resonance with the surrounding nonexcited 
extended lattice modes, and is also known as selftrap- 
ping. 

The normal modes (NMs) of a one-dimensional lin- 
ear chain with uncorrelated random potential are spa- 
tially localized (Anderson localization). Therefore any 
wavepacket, which is initially localized, remains localized 
for all time Anderson localization is harvesting on 
destructive interference and phase coherence, despite the 
fact, that the frequency of a localized NM is not tuned 
out of resonance with other NMs. From the experimen- 
tal point of view, disorder can be easily implemented in 
BEC's by tuning the wavelength of laser beams Q or 
by superimposing a speckle pattern [6] . Recently experi- 
ments have been reported, which observe Anderson local- 
ization for nonintcracting BECs In optics, disorder 
can be implemented in waveguide arrays by altering the 
fabrication process [1, Q or, in optically induced lattices, 
by adding a speckle beam [l0| . 

Nonlinearity induces interaction between NMs, and 
frequency shifts. While interaction favours delocaliza- 
tion [ill, [l2, [l^l I frequency shifts may lead to selftrap- 
ping and again to localization [l3| ■ Continuation of NMs 
into the nonlinear regime may keep localization, but also 
delocalize excitations via resonances flE]. Experimental 
studies of the combined action of nonlinearity and disor- 
der are possible both in BEC systems, as well as in optics 

According to Ref. [ij], any initial wavepacket and a 
fixed value of the nonlinearity parameter define three 
regimes - weak, intermediate and strong nonlinearity. 



The weak nonlinearity regime is characterized by An- 
derson localization on potentially large time scales, and 
subsequent detrapping and spreading. The intermediate 
regime is yielding spreading from scratch. The spreading 
continues despite the fact that the nonlinear frequency 
shifts weaken, since that is balanced by the increase in 
the number of excited NMs. Finally the strong nonlinear- 
ity regime leads to partial selftrapping, i.e. a part of the 
wavepacket selftraps and does not delocalize while the 
remaining part spreads again. The spreading is universal 
and subdiffusive, therefore rather slow, posing challenges 
for experimental studies, where one has to compete with 
dissipation mechanisms which lead to dephasing. 

Here we study the spreading of a wavepacket in a fi- 
nite lattice by ramping the strength of nonlinearity in 
time. The increase of nonlinearity with time counteracts 
the above diminishing of nonlinear frequency shifts, and 
substantially speeds up the spreading process. At the 
same time, resonant adiabatic excitation of distant NMs 
can contribute to a spreading as well. That makes our 
results also easier accessible for experiments. 



II. MODEL 

A. Equations of motion 

We consider a discrete nonlinear Schrodingcr (DNLS) 
model which describes the propagation of light in non- 
linear waveguide arrays or the evolution of a BEC in a 
periodic optical potential P, 0] . The Hamiltonian 
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where tpi{z) is a complex wave amplitude at site /. z 
corresponds to the dynamical variable (propagation co- 
ordinate for photons or time for atoms). The nonlinear 
cocfBcient P = (3{z) = cqz'^'- increases with z (/i > 0). We 
will numerically investigate linear ramping fi — 1 with 
velocity cq > 0. This can be implemented for the case of 
a BEC, where the interaction between condensed atoms 
can be described by a single parameter, the scattering 
length, which is proportional to an external applied mag- 
netic field. This parameter can be adjusted via Feshbach 
resonances and therefore the nonlinear interaction be- 
tween particles can be tuned and controlled in time ■ 
In optical systems, the nonlinearity can be adjusted in 
the propagation direction (z) in the fabrication process 
of laser-written waveguide arrays or by controlling 
the doping concentration in photovoltaic samples 19]. 

The random on-site energies e; are chosen uniformly 
from the interval [— ^] , where W is the strength of 
disorder. The norm A/" = is dynamically con- 

served, i.e. its value does not change with z. In our 
calculations Af — 1. We note that varying the norm 
is strictly equivalent to varying /3. For cq — 0, ^ 
is reduced to the linear eigenvalue problem ~Xi,A^j — 
+ + The eigenvectors A^,i are the 

NMs and the eigenvalues Ai^ their frequencies. 

We analyze normalized distributions ni — > us- 
ing the participation number P = 1/ nf, which mea- 
sures the number of the strongest excited sites, and the 
second moment TO2 = ~ which measures the 

size of the wavepacket. Here I = J2i ^"-i- 

B. The linear case co = 

Let us first consider an infinite system. For W = 
Cq ~ 0, solutions of Eq. IH) are extended plane waves 
ui{z) = uaexp[i{kl — Xkz)] with the dispersion relation 
Afe = —2 cos k. The frequencies are confined to the in- 
terval [—2, 2] {k is related to the input angle for optics or 
to the quasi-momentum for BEC). The group velocities 
of the plane waves \vg\ < 2. 

For W ^ 0,co = 0, the frequency spectrum is confined 
to [-2 - , 2 ]. The width of the spectrum {A^.} is 
A — A + W. All eigenvectors will be localized in space 
The asymptotic spatial decay of an eigenvector is given 
by A^^i - e-«' where ^(A^) < ^(0) « lOO/VF^ is the 
localization length The NM participation number 

Piy = ^/J2i^ti characterizes the spatial extend - local- 
ization volume - of the NM. It is distributed around the 
mean value p « 3.6^ with variance ~ (1-3^)^ [llj. The 
average spacing of eigenvalues of NMs within the range 
of a localization volume is therefore A A w A/p. 

The linear case is characterized by two frequency 
scales: the width of the spectrum A, and the average 
spacing AA < A. In addition, since we will operate with 
finite systems, we also have two spatial scales - the local- 
ization length ^ and the system size N. 

We launch a single site excitation in the center of our 



system with N = 100, and follow its evolution until 
z = Zmax {zmax = N / A) . For W = the fastest plane 
waves will reach the boundaries exactly at z — z^ax- 
A plot of the corresponding norm profile is shown in 
Figllja). The corresponding participation number is 
P K, AA — 0.44A^. Indeed, an almost completely spread 
wavepacket is characterized by roughly P = N /2, since a 
strictly homogeneous distribution is counterbalanced by 
dynamical fluctuations, and therefore, on average, every 
second site is not excited. Now we increase the strength 
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FIG. 1; Linear propagation. Spatial profile at z — Zmax 
for (a) W ^ 0, (h) W = 0.5 and (c) W = 3. N = 100. 

of disorder W. For W = 1 the localization length equals 
the system size. Therefore, for W = 0.5 [Figlljb)] we 
still observe almost full spreading, although the fronts 
propagate slower. For W = 3 [Figdjc)] we observe An- 
derson localization - the wavepacket is confined to the 
localization volume which is of the order of 20 sites. 

We performed runs for 100 realizations, and compute 
the average of the second moment m2 and of the par- 
ticipation number P at z^ax (see Figl2]). The second 
moment m2 decreases with increasing disorder as ex- 
pected, showing that waves propagate slower even within 
the localization volume {W < 1) due to disorder induced 
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backscattering. However, the participation number P 
shows a shght increase for weak disorder, with a subse- 
quent (expected) decrease for stronger disorder. For the 
ordered case W — the chosen initial condition excites 
one half of all available NMs. That follows from the re- 
flection invariance of the lattice around the center, where 
the initial wavepacket is placed. All NMs separate into 
two irreducible groups of even and odd ones, or symmet- 
ric and antisymmetric ones with respect to reflections 
around the chain center. A single site excitation in the 
center excites only even NMs. For weak disorder all NMs 
start to become excited. Therefore the part of the vol- 
ume which is occupied by the wavepacket at some later 
time, is filled more homogeneously. Note that the effect 
is weak - the participation number increases by ~ 15%. 
That may be due to the fact, that at the same time the 
wavepacket spreads less effectively, as it follows from the 
results on the second moment. 
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FIG. 2: Linear propagation, (a) Average second moment 
m2 measured at time Zmax as a function of W. (b) Same 
as (a), but for the participation number P. N = 100 and 
averaging was performed with 100 realizations. 



C. The nonlinear case 

The equations of motion in normal mode space read 
- i^^ ^ X^cf)^ + f3 ^ /i/, 1^1,1/2, '/'i^iCa'^f 3 (2) 

with the overlap integral 

I 



The variables (^y determine the complex time-dependent 
amplitudes of the NMs. 

Nonlinearity therefore induces interaction between 
NMs. In particular the nonlinearity renormalizes fre- 
quencies. If frequencies are shifted out of the spectrum of 
the linear equation, selftrapping occurs, and excitations 
stay localized for long, may be infinite, times. We are not 
aware of a straightforward and unique way to calculate 
such a frequency shift for a given distribution n; . There- 
fore we will look for suitable estimates. One possibility is 
to neglect the coupling along the lattice, and treat lattice 
sites as independent. Then the nonlinear frequency shift 
at site Hs (5A « —/3ni. Another possibility is to derive 
an effective frequency Ae// for a given state. We treat a 
state as a stationary one i^ii^) = ^/ exp (— ?Ae//z). We 
insert this ansatz in ([1]), multiply this equation by ipi 
and sum over all lattice sites: 

^e// = - E h I' + V'^+i V'r + i^t+ii^i (4) 
I 

At any time z the effective frequency Xeff{z) = Hiz) — 
l3{z)/[2P{z)]. In practice both definitions from above are 
giving similar results, and we will mainly use 

In all our simulations we use a single-site excitation as 
the initial condition: tpi{0) — Si,i^ with Ic — N/2, such 
that A(0) = -e/^ and P(0) = 1. To effectively harvest 
on resonances with the spectrum we have imposed that 
e;e — —W/2. With that, and for cq > 0, the effective 
frequency will decrease starting from Ae// = W/2 until 
it reaches the bottom of the spectrum X^ot where we stop 
our simulations. The bottom of the spectrum of an infi- 
nite system is located at A = —2 — W/2. However, for fi- 
nite systems the bottom of the spectrum Xbot > —2— W/2 
depends on the given realization. From ([4]) we see that 
a positive increment in P will decrease the effective fre- 
quency of the system. Therefore, the state will be able to 
resonate with other NMs inside the spectrum. Outside 
of it, the solution transforms into a selftr app ed localized 
state similar to a discrete soliton [1, which is a 

time-periodic and exponentially localized excitation p|. 

III. NUMERICAL RESULTS 

In Fig. [3] we show the evolution of an initial single site 
excitation for W = 3 and N = 100, when < N. For 
Co = 0, the excitation remains trapped due to Anderson 
localization, and does not reach the boundaries of the fi- 
nite chain [see Fig. [2l^c)]. However, for co = 0.2, a slow 
increase of nonlinearity leads to a complete delocalization 
of the wavepacket at 2; « 800. Indeed, the effective fre- 
quency A decreases and around that time reaches the bot- 
tom of the spectrum. At the same time, the participation 
number reaches its saturation value P « N/2 — 50. For 
a larger ramping velocity cq = 1, but exactly the same 
disorder realization, the wavepacket does not spread over 
the entire chain. At a much shorter time z w 70 the effec- 
tive frequency touches the bottom of the spectrum, the 
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participation number saturates around P « 20, and the 
state becomes selftrapped, still occupying only 20 out of 
100 sites. Further increase of the ramping velocity will 
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FIG. 3: Nonlinear propagation, (a) X^ff versus /3; (b) P 
versus 13; (c) final profiles. Gray and black curves correspond 
to CO = 1 and 0.2, respectively. W = 3 and iV = 100. The 
horizontal dashed line in (a) corresponds to Xtot- 

make the final wavepacket more and more localized. 

In order to study the dependence of the delocaliza- 
tion process on the disorder strength, we define our sim- 
ulation scheme as follows: (i) we ramp the nonlinear- 
ity for a given velocity cq; (ii) we stop the simulation 
when Xeff reaches the bottom of the spectrum Xbot, in- 
dicating a stop of spreading due to selftrapping; (iii) we 
compute P and if P < 40 we decrease the velocity cq 
until P > 40, which corresponds to a fully delocalized 
wavepacket. Therefore we obtain the critical velocity 
Ccr for a given disorder realization. For cq > Ccr the 
wavepacket does not spread over the entire chain, while 
for Co < Ccr it does. We repeat the scheme for 100 differ- 
ent disorder realizations and obtain an average value for 
c^. Finally we change the disorder strength W and ob- 
tain the dependence c^{W). Results are shown in FigH) 




FIG. 4: Critical velocity versus degree of disorder for 100 
realizations. Line represents the average tw. A'' = 100. 

We observe two different regimes. For strong disorder 
W > 1 the localization length is smaller than the system 
size N — 100. The critical velocity is monotonously 
decreasing with increasing disorder strength. This can 
be expected, since the stronger the disorder, the more lo- 
calized the NMs are, and the less is the number of other 
NMs a given NM can interact with due to nonlinearity. 
Therefore we need slower ramping in order to give ample 
time to NMs to interact. For weak disorder W < 1 the lo- 
calization length is larger than the system size N — 100. 
Therefore all NMs are essentially extended over the en- 
tire chain. In this regime, we observe that the critical 
velocity increases with increasing disorder strength. 
This effect may be due to the fact, that in the limit of 
zero disorder W = 0, the nonlinearity induces a selective 
interaction between NMs, as mentioned above (for de- 
tails of that selective interaction, see Ref . [2^ ) . Weak but 
nonozero disorder removes the selection rules, and leads 
to an interaction of all modes with each other, thereby 
increasing the number of available modes (channels) by a 
factor of two. The nonlinearity induced spreading is more 
effective. In particular it can redistribute the energy into 
two times more modes. 

The hallmark of the transition from weak to strong 
disorder is the maximum in the dependence of c^(VF) at 
W = 1. The maximum is located at W ^ ^/^00/N. In 
particular, for = 50 it is located at W ^ 1.4 and for 
iV = 200 at PF 0.7. For large TV > ^ the maximum 
position is therefore shifted closer to VF = 0, and the 
value in the maximum increases. The value of at 
W = can be estimated by noting that the fastest plane 
waves reach the edges of the chain at a time zi, w iV/4. 
The time it takes to shift the frequency of the initially 
excited oscillator out of the band is Zg ~ 4/co. Equating 
Zh = Zs we estimate c^{W = 0) ~ 16/A^. The increase 
of this critical speed in the presence of weak disorder is 
due to the above discussed increase of modes (channels) 
available. This increase by a factor of two leads to a 
decrease of the norm in each mode also by a factor of 
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two, and, consequently, implying a larger critical velocity. 
Therefore ^ 2/co and max(c^(M^ 7^ 0)) ~ 2(^{W = 
0) at best. This increase by a factor of two is close to the 
numerical observation (see FigH]). 



IV. DIFFUSION VERSUS RESONANT 
SPREADING 

Let us discuss possible mechanisms of wavepacket 
spreading in the regime of strong disorder, when < N. 
The packet spreading beyond the localization volume is 
entirely due to the nonlinearity, which induces interac- 
tions between NMs. Assume the wavepacket has a certain 
size at some time z. Further spreading implies excitation 
of exterior NMs. Since the interaction between NMs falls 
off exponentially for distances larger than ^, the relevant 
exterior NMs to be excited will be located in a nonex- 
cited (cold) boundary layer outside the wavepacket, with 
a width roughly of the order of the localization volume. A 
given cold exterior NM can be excited coherently or inco- 
herently. A coherent excitation implies a resonance with 
an excited NM from the wavepacket, and a correspond- 
ing resonant transfer of energy (as it happens during the 
beating of energy between two weakly coupled harmonic 
oscillators) . An incoherent excitation implies the absence 
of such resonances, and an almost random fluctuation of 
the NMs phases inside the wavepacket. It generates a 
random force which incoherently excites (heats) the ex- 
terior NMs. Such a spreading corresponds to a diffusive 
spreading of the wavepacket, therefore of the norm, and 
the energy. 



A. Diffusion 

The spreading of wavepackets in the presence of a con- 
stant nonzero nonlinearity strength (3 = Pq was studied 
in several papers [HI, [H, [11, [13] ■ In particular, it was ob- 
served, that a spreading wavepacket is characterized by 
a growth of its second moment 7712 ~ z". The exponent 
a = 0.33 ± 0.02 [l^. A theoretical analysis showed, that 
this subdiffusive spreading is not due to resonant exci- 
tation of exterior modes, but due to incoherent heating. 
The origin of the chaotic dynamics of the wavepacket it- 
self comes from resonant interaction of NM pairs inside 
the packet. The statistical analysis of the probabilities 
of such internal resonances leads to the conclusion, that 
the second moment obeys the following equation: 



dm 2 



(5) 



Here n is the average norm density inside the wavepacket, 
and C{W) an unknown function, which however, as nu- 
merical studies suggest, decreases with increasing W. 
Since the packet size is ^ it follows that the par- 
ticipation number P ^ 1/n and the second moment 



1712 l/n"^. As a consequence, m2 ~ {P^zY^^. The ex- 
ponent a = 1/3 and is in very good agreement with the 
numerical studies. The nonlinearity induced frequency 
shift inside the packet S\ ~ /3n ^ z~^l^ . The more 
the packet spreads, the smaller the frequency shifts are. 
This weakening is counterbalanced by the increase in the 
size of the packet, so that more modes are involved, and 
guarantee a slow but steady subdiffusive spreading. 

Assuming, that in the present case the spreading is 
incoherent as well, we find with /3 = cq^^: 



m2 



4/3 (4p+l)/3 



2/3 (4^+l)/6 



The nonlinearity induced frequency shift 



(6) 



(7) 



In the present numerical studies /i = 1 , and we find 7712 ~ 
2;5/3^ p z^l^ and (5A ~ z^l^ . At variance with the 
case of constant nonlinearity (which yields subdiffusion 
with a = 1/3) we obtain superdiffusion with a — 5/3. 
This is due to the fact, that the increasing nonlinearity 
counterbalances the decay of the frequency shifts, and 
enhances the interaction of the NMs from the wavepacket 
with cold exterior NMs. The growth of the frequency 
shifts will finally lead to a selftrapped state, as observed 
in the numerical studies. 

Let us estimate the dependence of the participation 
number P of the selftrapped wavepacket on cq- With 
Eq.® it follows 



-2/3 (4p+l)/6 



Therefore the nonlinear frequency shift &\ 



(8) 



Zq^^z^"^^ When this frequency shift reaches some fi- 

- Zf, the wavepacket 

Finally 



nite (disorder dependent) value at z 
selftraps. Therefore Zf 



-2/(2m-1) 



1(2/) 



^-1/(2^-1) 
-0 



(9) 



For /i = 1 we find that the selftrapped wavepacket has 
size P ^ I/cq. Increasing the ramping speed cq therefore 
leads to a smaller extend of the selftrapped wavepacket, 
as observed in the numerical results. We can not estimate 
the measured dependence c^(VF), since we do not know 
the function C{W) from (P. 



B. Resonant spreading 

Let us assume that the wavepacket spreads by reso- 
nantly exciting exterior NMs. Except for exponentially 
small ramping velocities, the relevant exterior NMs will 
be located in a layer of the width of the localization vol- 
ume. The average frequency spacing A A « A/p sepa- 
rates possible frequencies of resonant interactions, which 
can take place due to the nonlinearity induced frequency 
shift. Therefore the number of possible resonances is 
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limited to the number of NMs within one locahzation 
volume, i.e. to the localization volume p itself. The 
maximum size, to which a wavepacket can grow, is then 
proportional to the p^. The numerical results in Fig[3] 
show, that for strong disorder W — 7 the critical veloc- 
ities are of the order of w 0.01. For such a velocity 
the wavepacket spreads completely over a chain with 100 
sites. The maximum localization length ^ ~ 2, therefore 
the localization volume is less than p « 7, and the max- 
imum size to less than p^ « 50, which indicates, that 
resonant spreading is not the dominant mechanism. 

Let us assume the optimum case scenario for resonant 
spreading. That implies, that at a given time a resonance 
at the edge of the wavepacket can take place, and the 
energy is transfered much faster, than it takes to shift the 
frequency to the next resonance. Assume we start with 
one NM excited, with its frequency located at the edge 
of the spectrum. For /i = 1, we will hit the next possible 
resonance after time zi — AX/cq. We assume, that a 
resonant mode is available in the localization volume, 
and is immediately excited. The norm is now equally 
distributed between both NMs, each of them having norm 
1/2. Nonlinearity is further ramped up, but in order to 
shift the frequency to the next resonance, a larger time 
Z2 — 2zi is needed, and so on. As long as the wavepacket 
does not selftrap, we find Zj = jzi. Therefore the total 
time to reach the j-th resonance scales as z ~ j^, the 
participation number P ^ j grows as P ^ \/z and the 
second moment m2 ~ z. Since that idealized process 
is already slower than the incoherent spreading due to 
diffusion, we expect that resonant spreading is weakly 
contributing to the numerically observed spreading. 

V. CONCLUSIONS 

We have investigated the spreading of a wavepacket in 
a disordered nonlinear chain, when ramping the nonlin- 



earity in time. For linear ramping /i = 1 we find that the 
spreading is mostly due to incoherent diffusive processes. 
For an infinite chain, the nonlinear frequency shift always 
leads to a selftrapping of the wavepacket, and therefore 
not to a complete delocalization. The second moment 
1712 ^ z^/"^ is predicted to grow superdiffusively, at vari- 
ance to the previously studied case of constant nonzero 
nonlinearity, which yields subdiffusion. Therefore the 
present case is easier accessible in experiments, due to 
possible restrictions on maximum propagation times. Ex- 
periments are done with finite systems. We studied the 
case of a finite chain, and computed the critical ramp- 
ing speed Ccr, which separates selftrapped localized from 
selftrapped delocalized states at the end of the ramping 
process. We find that Ccr grows with increasing strength 
of disorder when the localization length is larger than 
the system size. That increase is due to the fact, that 
disorder is removing selection rules for the interaction of 
normal modes, present in the case of zero disorder. For 
strong disorder, i.e. when the localization length becomes 
smaller than the system size, the critical velocity drops 
with further increase of the disorder strength. That is 
due to the fact, that the normal modes are less and less 
spread, interact with fewer other modes, and the nonlin- 
earity induced frequency shift is more efficient in tuning 
the wavepacket faster into a selftrapped state. 
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